home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / prog / yamp2.zip / TESTGRAF.CPP < prev    next >
C/C++ Source or Header  |  1992-01-16  |  5KB  |  193 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/23/92
  7.  
  8. Copyright(c) Mark Von Tress 1992
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20. #include "virt.h"
  21.  
  22. //required global declaration for the
  23. //  matrix stack object
  24.  
  25. // unsigned int _stklen = STACKLENGTH;
  26.  
  27. MStack *Dispatch = new MStack;
  28.  
  29.  
  30. VMatrix &getx( int N )
  31.   // create an x matrix
  32.   {
  33.       Dispatch->Inclevel();
  34.       VMatrix x, c1, x2;
  35.  
  36.       c1 = Fill(N,1,1.0);
  37.       x = Index( 1, N ) - ((double)N)*0.5;
  38.       x = Ch( c1,x );
  39.  
  40.       // push x onto the stack
  41.       Dispatch->Push(x);
  42.       // decrement the subroutine nesting level
  43.       // and return the stack top
  44.       return Dispatch->DecReturn();
  45.   }
  46.  
  47. VMatrix &gety( VMatrix &x, VMatrix &beta)
  48.   // create a y matrix
  49.   {
  50.       Dispatch->Inclevel();
  51.       VMatrix y;
  52.  
  53.       y = x*beta;
  54.       srand(123);
  55.       for(int i=1; i<=y.r; i++) {
  56.         // use sum of 3 uniforms for an approximate
  57.         // normal random variable
  58.         int u = random(100)+random(100)+random(100)+3;
  59.         y.M(i,1) = y.m(i,1) + 5.0*sin( 3.14*((double) (i%4))/3.0 )
  60.                  + ((double) (u-150))/300.0;
  61.       }
  62.       Dispatch->Push(y);
  63.       return Dispatch->DecReturn();
  64.   }
  65.  
  66. VMatrix ®ression( VMatrix& x, VMatrix& y )
  67. // do a multiple linear regression
  68.       {
  69.       Dispatch->Inclevel();
  70.       VMatrix yx, reg, betahat;
  71.       int N=x.r, p=x.c;
  72.  
  73.       // solve for regression parameters using sweep
  74.       yx = Ch(y,x);
  75.       reg = Sweep( 2,p+1, Tran(yx)*yx);
  76.       // calculate mean square error
  77.       reg.M(1,1) = reg.m(1,1)/((double )( N-p));
  78.       reg.DisplayMat();
  79.  
  80.  
  81.       // solve regression using normal equations
  82.       betahat = Inv(Tran(x)*x)*Tran(x)*y;
  83.  
  84.       Dispatch->Push( betahat );
  85.       return Dispatch->DecReturn();
  86. }
  87. void plotResiduals( VMatrix &resids )
  88.   {
  89.     Dispatch->Inclevel();
  90.     VMatrix grf = Ch( Index( 1, resids.r ), resids );
  91.     GMatrix Agraph( grf, -'%' );
  92.     *Agraph.PathToDriver = "C:\\tc\\bgi";
  93.     *Agraph.title = "Residuals for data";
  94.     *Agraph.title2= "with serial correlations with frequency 0.25";
  95.     *Agraph.yname = "Residuals";
  96.     *Agraph.xname = "Index";
  97.     Agraph.Href( 0.0 );
  98.     Agraph.Show();
  99.     Dispatch->Cleanstack();
  100.     Dispatch->Declevel();
  101.   }
  102.  
  103.  
  104. VMatrix &GetSerialCovar( VMatrix &R, VMatrix &spectdens )
  105.   {
  106.     Dispatch->Inclevel();
  107.     VMatrix centered, z, covar;
  108.     double n = (double) R.r;
  109.     centered = R - Sum( R/n ).m(1,1);
  110.     z = Fft( centered );
  111.     spectdens = Sum( z%z/n, COLUMNS);
  112.     covar = Fft( spectdens, FFALSE );
  113.     Dispatch->Push( covar );
  114.     return Dispatch->DecReturn();
  115.   }
  116. void plotCorrelogram( VMatrix &serial )
  117.   {
  118.      Dispatch->Inclevel();
  119.  
  120.      int n = serial.r;
  121.      double sigma = serial.m(1,1);
  122.      VMatrix Correlogram = serial/sigma;
  123.      Correlogram = Submat( Correlogram, n/2, 1, 2, 1);
  124.      VMatrix graf = Ch( Index( 1, Correlogram.r ), Correlogram );
  125.  
  126.      GMatrix Agraph( graf );
  127.      double dn2 = 2.0/sqrt((double)n);
  128.      *Agraph.PathToDriver = "C:\tc\bgi";
  129.      *Agraph.title = "Serial Correlations";
  130.      *Agraph.title2= "for sample residuals";
  131.      *Agraph.yname = "Serial correlations";
  132.      *Agraph.xname = "Lags";
  133.      Agraph.Href( 0.0 );
  134.      Agraph.Show();
  135.  
  136.      Dispatch->Cleanstack();
  137.      Dispatch->Declevel();
  138.   }
  139.  
  140. void plotPeriodogram( VMatrix &spectdens )
  141.   {
  142.      Dispatch->Inclevel();
  143.      int n = spectdens.r;
  144.      double dn = (double)n;
  145.      double sigma = (Sum(spectdens).m(1,1) - spectdens.m(1,1));
  146.      VMatrix Periodogram = spectdens/sigma;
  147.      Periodogram = Mlog(Submat( Periodogram, 1+n/2, 1, 2, 1)) ;
  148.      VMatrix graf = Ch( Index( 1, Periodogram.r )/dn, Periodogram );
  149.  
  150.      GMatrix Agraph( graf );
  151.      *Agraph.PathToDriver = "C:\tc\bgi";
  152.      *Agraph.title = "Periodogram";
  153.      *Agraph.title2= "for sample residuals";
  154.      *Agraph.yname = "Log periodogram";
  155.      *Agraph.xname = "Frequencies";
  156.      Agraph.Vref( 0.25 );
  157.      Agraph.Show();
  158.  
  159.      Dispatch->Cleanstack();
  160.      Dispatch->Declevel();
  161.   }
  162.  
  163.  
  164. main()
  165. {
  166.       Dispatch->Inclevel();
  167.       VMatrix x, beta("beta",2,1), y, betahat, resids, serial;
  168.       Setwid(15);
  169.       Setdec(10);
  170.  
  171.       beta.M(1,1) = 1;
  172.       beta.M(2,1) = 0.5;
  173.  
  174.       x = getx( 128 );
  175.       y = gety(x,beta);
  176.  
  177.       betahat = regression(x,y);
  178.       betahat.Nameit( "Text book betahat");
  179.       (Tran(beta)).DisplayMat();
  180.       (Tran(betahat)).DisplayMat();
  181.  
  182.       resids = y - x*betahat;
  183.       plotResiduals( resids );
  184.  
  185.       VMatrix spectdens;
  186.       serial = GetSerialCovar( resids, spectdens );
  187.  
  188.       plotCorrelogram( serial );
  189.       plotPeriodogram( spectdens);
  190.  
  191.    vclose();
  192. }
  193.